clear; clc; close all;

%% 导入数据
load result05_10_with_IC_filtered.mat;
data = IC_data;

% 初始化
num_data = length(data);

% 用于存储特征值的数组
peak_features = zeros(num_data, 8); % 每个循环有8个特征值

for idx = 1:num_data
    
    % 选择要分析的数据
    x = data{idx};
    
    % 检查 x 的维度是否正确
    if isempty(x) || size(x, 2) < 2
        % 如果 x 为空或列数不足，跳过当前循环
        warning('数据编号 %d 的维度不正确，跳过。', idx);
        peak_features(idx, :) = NaN(1, 8); % 用 NaN 填充
        continue;
    end
    
    V = x(:, 1); % 第一列是电压
    dQdV_filter = x(:, 2); % 第二列是 dQ/dV 数据

    % 删除电压大于 3.9 的部分
    valid_indices = V <= 3.9; % 找到电压小于等于 3.9 的索引
    V = V(valid_indices); % 保留有效电压
    dQdV_filter = dQdV_filter(valid_indices); % 保留对应的 dQ/dV 数据

    %% 提取峰值特征
    [pks, locs, widths, prominences] = findpeaks(dQdV_filter, V,  ...
        'MinPeakHeight', 0.5, ...       % 最小高度'MinPeakProminence', 0.02, ...  % 最小突出度
        'MinPeakDistance', 0.1);        % 最小距离

    % 按电压位置对峰值排序（从小到大）
    if ~isempty(pks)
        [sorted_locs, sorted_idx] = sort(locs, 'ascend'); % 按电压位置排序

        % 提取峰值1（电压较小的第一个峰）
        if length(sorted_locs) >= 1
            peak1 = pks(sorted_idx(1)); % 峰值1大小
            loc1 = sorted_locs(1); % 峰值1位置

            % 计算峰值1的宽度
            left_index = find(dQdV_filter >= pks(sorted_idx(1)) / 2, 1, 'first');
            right_index = find(dQdV_filter >= pks(sorted_idx(1)) / 2, 1, 'last');
            peak_width1 = V(right_index) - V(left_index); % 峰值1宽度

            % 计算峰值1面积（使用梯形法则）
            valid_index1 = (V >= loc1 - peak_width1/2 & V <= loc1 + peak_width1/2); % 有效区间
            if any(valid_index1)
                x_valid1 = V(valid_index1);
                y_valid1 = dQdV_filter(valid_index1);
                main_peak_area1 = sum((y_valid1(1:end-1) + y_valid1(2:end)) / 2 .* diff(x_valid1)); % 梯形法则
            else
                main_peak_area1 = 0; % 如果没有有效区间，面积为0
            end
        else
            % 如果没有峰值1，相关数据设为0
            peak1 = 0;
            loc1 = 0;
            peak_width1 = 0;
            main_peak_area1 = 0;
        end

        % 提取峰值2（电压大于 3.5V 的峰）
        if length(sorted_locs) >= 2
            % 找到电压大于 3.5V 的峰
            peak2_candidates = sorted_locs(sorted_locs > 3.5);
            if ~isempty(peak2_candidates)
                % 选择电压大于 3.5V 的第一个峰
                peak2_loc = peak2_candidates(1);
                peak2_idx = find(sorted_locs == peak2_loc, 1);
                peak2 = pks(sorted_idx(peak2_idx)); % 峰值2大小
                loc2 = sorted_locs(peak2_idx); % 峰值2位置

                % 计算峰值2的宽度
                left_index = find(dQdV_filter >= pks(sorted_idx(peak2_idx)) / 2, 1, 'first');
                right_index = find(dQdV_filter >= pks(sorted_idx(peak2_idx)) / 2, 1, 'last');
                peak_width2 = V(right_index) - V(left_index); % 峰值2宽度

                % 计算峰值2面积（使用梯形法则）
                valid_index2 = (V >= loc2 - peak_width2/2 & V <= loc2 + peak_width2/2); % 有效区间
                if any(valid_index2)
                    x_valid2 = V(valid_index2);
                    y_valid2 = dQdV_filter(valid_index2);
                    main_peak_area2 = sum((y_valid2(1:end-1) + y_valid2(2:end)) / 2 .* diff(x_valid2)); % 梯形法则
                else
                    main_peak_area2 = 0; % 如果没有有效区间，面积为0
                end
            else
                % 如果没有电压大于 3.5V 的峰，相关数据设为0
                peak2 = 0;
                loc2 = 0;
                peak_width2 = 0;
                main_peak_area2 = 0;
            end
        else
            % 如果没有峰值2，相关数据设为0
            peak2 = 0;
            loc2 = 0;
            peak_width2 = 0;
            main_peak_area2 = 0;
        end
    end
    
    % 存储特征值（按指定顺序）
    peak_features(idx, :) = [peak1, loc1, peak_width1, main_peak_area1, peak2, loc2, peak_width2, main_peak_area2];
end


% 对噪声处理后的数据进行高斯滤波
w_gauss = gausswin(50, 3); % 高斯滤波参数
w_gauss = w_gauss / sum(w_gauss);
peak_features_filtered = zeros(size(peak_features));
for col = 1:8
    % 获取当前列的特征值
    current_feature = peak_features(:, col);
    
    % 检查是否存在非有限值（NaN 或 Inf）
    if any(~isfinite(current_feature))
        % 如果有非有限值，跳过滤波并保留原始值
        warning('列 %d 包含非有限值，跳过滤波。', col);
        peak_features_filtered(:, col) = current_feature;
    else
        % 如果数据是有限的，进行高斯滤波
        peak_features_filtered(:, col) = filtfilt(w_gauss, 1, current_feature);
    end
end

% 保存噪声处理前、噪声处理后和滤波后的数据到同一个 .mat 文件
save('peak_features05_10.mat', 'peak_features', 'peak_features_filtered');

% 画图表现随着数值编号的增加峰特征的变化趋势
figure;
hold on;
plot(peak_features_filtered(:, 1), 'DisplayName', '峰值1大小');
plot(peak_features_filtered(:, 2), 'DisplayName', '峰值1位置');
plot(peak_features_filtered(:, 3), 'DisplayName', '峰值1宽度');
plot(peak_features_filtered(:, 4), 'DisplayName', '峰值1面积');
plot(peak_features_filtered(:, 5), 'DisplayName', '峰值2大小');
plot(peak_features_filtered(:, 6), 'DisplayName', '峰值2位置');
plot(peak_features_filtered(:, 7), 'DisplayName', '峰值2宽度');
plot(peak_features_filtered(:, 8), 'DisplayName', '峰值2面积');
legend();
xlabel('数据编号');
ylabel('特征值');
title('峰特征随数据编号的变化趋势（电压 ≤ 3.9）');
hold off;